%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot_aggregate_vix.m
% plot the aggregate quantities and VIX 
%
% Ivan Alfaro, Nick Bloom and Xiaoji Lin 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all; clc; clearvars;

%Here you pre-define the sample that you want to use
year0=1991; quarter0=year0*10+1; month0=year0*100+01; %Should start in first quarter/month.
year1=2019; quarter1=year1*10+4; month1=year1*100+12; %Should end in last quarter/month.

%A description of data sets: missing value is set to zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data1: NIPA1.1.4: contains price indexes
% Data2: NIPAS1q: contains NIPA Integrated Macroeconomic Accounts for Total Economy
% Data3: NIPAS5q: contains NIPA Integrated Macroeconomic Accounts for Nonfinancial Corporate Business 
% Data4: vol: contains monthly VIS. 
% Data5: ShillerMonthly: contains monthly dividend from Shiller
% Data6: FlowFundF102q: contains Flow of Funds talbe 102 quarterly
% Data7: BLSEmploymentQ: contains quarterly private employment from BLS
% Data8: NIPATable1.1.5Annual.txt:          contains annual nominal annual private fixed nonresidential investment
% Data9: NIPATable1.1.5Quarterly.txt:       contains quarterly nominal annual private fixed nonresidential investment
% Data10: NIPATable1.1.9Quarterly.txt:       contains quarterly implicit price delfator of nonresidential investment
% Data11: NIPATableFixedAsset1.1Annual.txt:  contains annual private nonresidential capital stock at year-end prices
% Data12: NIPATableFixedAsset1.3Annual.txt:  contains annual nominal annual private nonresidential depreciation
% Data13: NIPATable1.1.9Annual.txt:          contains annual implicit price delfator of nonresidential investment

% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Data1 = load('Data\NIPA1.1.4q_new.txt');                            ini1=find(Data1(:,1)==quarter0);    fin1=find(Data1(:,1)==quarter1); 
Data2 = load('Data\NIPAS1q_new.txt');                               ini2=find(Data2(:,1)==quarter0);    fin2=find(Data2(:,1)==quarter1);                 
Data3 = load('Data\NIPAS5q_new.txt');                               ini3=find(Data3(:,1)==quarter0);    fin3=find(Data3(:,1)==quarter1);
Data4 = load('Data\Vol_monthly.txt');                               ini4=find(Data4(:,1)==month0);      fin4=find(Data4(:,1)==month1);
Data5 = load('Data\ShillerMonthly_new.txt');                        ini5=find(Data5(:,1)==month0);      fin5=find(Data5(:,1)==month1);
%Data6 = load('Data\FlowFundF102q_new.txt');                        ini6=find(Data6(:,1)==quarter0);    fin6=find(Data6(:,1)==quarter1); 
Data7 = load('Data\BLSEmployment_monthly.txt');                     ini7=find(Data7(:,1)==month0);      fin7=find(Data7(:,1)==month1); 
Data8 = load('Data\NIPA1.1.5a_new.txt');                            ini8=find(Data8(:,1)==year0);       fin8=find(Data8(:,1)==year1); 
Data9 = load('Data\NIPA1.1.5q_new.txt');                            ini9=find(Data9(:,1)==quarter0);    fin9=find(Data9(:,1)==quarter1);
Data10 = load('Data\NIPA1.1.9q_new.txt');                           ini10=find(Data10(:,1)==quarter0);  fin10=find(Data10(:,1)==quarter1);
Data11 = load('Data\NIPATableFixedAsset1.1a_new.txt');              ini11=find(Data11(:,1)==year0);     fin11=find(Data11(:,1)==year1);    
Data12 = load('Data\NIPATableFixedAsset1.3a_new.txt');              ini12=find(Data12(:,1)==year0);     fin12=find(Data12(:,1)==year1);    
Data13 = load('Data\NIPA1.1.9a_new.txt');                           ini13=find(Data13(:,1)==year0);     fin13=find(Data13(:,1)==year1);    

% calculate inflation
Cpi   = Data1(:,3);                             % CPI index

% GDP
dr      = 100;
GDP     = Data2(:,2);
GDPR    = GDP(ini2:fin2)./Cpi(ini1:fin1);
GDPGr   = GDP(ini2:fin2)./GDP(ini2-1:fin2-1)./((Cpi(ini1:fin1)./Cpi(ini1-1:fin1-1))) - 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% flow of funds data
DebtST = Data3(:,139)+Data3(:,143);  % open market paper + short-term loans = short-term debt
DebtLT = Data3(:,141)+ Data3(:,147); % bonds + mortgages = long-term debt
Debt   = DebtST+DebtLT;              % total debt
DebtGr   = Debt(ini3:fin3)./Debt(ini3-1:fin3-1)./((Cpi(ini1:fin1)./Cpi(ini1-1:fin1-1))) - 1;

% liquid asset or flow
LiquidAT = Data3(:,110);              % liquid assets: Time and savings deposits 
LiquidATGr = LiquidAT(ini3:fin3)./LiquidAT(ini3-1:fin3-1)./((Cpi(ini1:fin1)./Cpi(ini1-1:fin1-1))) - 1;

% real dividend from Shiller
Divm  = Data5(:,2);
Divq  = monthlyToQuarterly(Divm(ini5-3:fin5));
DivGr = [Divq(2:end)./Divq(1:end-1)]-1;

% employment from BLS
Empm   = Data7(:,4);
Empq  = monthlyToQuarterly(Empm(ini7-3:fin7));
EmpGr = Empq(2:end)./Empq(1:end-1)- 1;

%quarterly investment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Year    = Data8(ini8:fin8,1);
% Quarter = Data9(ini9:fin9,1);
Ia      = Data8(ini8:fin8,10);                                  % nominal annual I
Ka      = Data11(ini11:fin11,5);                               % nominal annual K
Da      = Data12(ini12:fin12,5);                               % nominal annual depreciation
Iq      = Data9(ini9:fin9,10);                                 % nominal quarterly I
Pq      = Data10(ini10:fin10,10);                              % quarterly implicit price deflator of investment
Pa      = Data13(ini13:fin13,10);                              % annual implicit price deflator of investment

% get quarterly adjusted investment following BCE
Iqq     = quarterlyToQuarterly(Iq,Ia);                      % following BCE
Iqq_r   = Iqq./Pq;                                          % real quarterly investment
 
% Inflation    
Inf_temp        = zeros(size(Data10,1),1);
Inf_temp(2:end) = Data10(2:end,10)./Data10(1:end-1,10)-1;             % quarterly inflation
Inf             = Inf_temp(ini10:fin10);

% get quarterly depreciation
T   = size(Pq,1);
k   = T/4;
Dq  = getQuarterlyDep(Da,Inf,k);

% to get quarterly capital
Kq    = zeros(T,1);
Kq0   = 2*Data10(ini10-1,10)*Data11(ini11-1,5)/(Data10(ini10-1,10)+Data10(ini10,10));          % the initial value of quarterly capital is the previous year-end annual capital
% Kq0   = Data4(ini4-1,5);          % the initial value of quarterly capital is the previous year-end annual capital
Pend  = Data10(fin10+1,10);
Kq    = getQuarterlyCap(Ka,Inf,Dq,Iqq, Pq, Kq0, Pend,T);

% Real quaterly capital and IK
Kq_r        = Kq./Pq;
IKq         = zeros(T,1);
IKq(2:end)  = Iqq_r(2:end)./Kq_r(1:end-1);
IKq(1)      = mean(IKq(2:4));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% volatility
Volm  = Data4(:,4);  
Volq  = monthlyToQuarterly(Volm(ini4:fin4));

%% figure
if 1
figure(2)

MktVol = (Volq)./100 - mean((Volq)./100);

date2 = year0:0.25:year1+1;
date  = date2(2:end-1);
fs   = 12;
rolN = 4;
ms   = 2;
time0 = 1990;
time1 = 2020;

subplot(3,2,1)
[AX,H1,H2] = plotyy(date,MA(IKq(2:end),rolN,1)-mean(MA(IKq(2:end),rolN,1)) ,date, MktVol(1:end-1),'plot');
set(H1,'LineStyle','-','LineWidth',ms,'color','b')
set(H2,'LineStyle','-.','LineWidth',ms,'color','r')
set(AX(1),'XLim',[time0 time1])
set(AX(2),'XLim',[time0 time1])
xmin = -.01;
xmax = 0.01; %round(max(MA(IKq(2:end),rolN,1)),2);
ymin = -.5;
ymax = 0.5;% round(max(MktVol(1:end-1)),2);
set(AX(1),'YLim',[xmin xmax])
set(AX(2),'YLim',[ymin ymax])
set(AX(1),'YTick',[xmin 0.5*(xmin+xmax) xmax]); 
set(AX(2),'YTick',[ymin 0.5*(ymin+ymax) ymax]); 
set(AX(1),'ycolor','b','Fontsize',fs) 
set(AX(2),'ycolor','r','Fontsize',fs)
lgd=legend( 'Real and financial variables','VIX implied S&P500 volatility');
set(lgd,'FontSize',fs)
legend('boxoff')
title('Investment rate','FontSize',fs) 
grid on


subplot(3,2,2)
[AX,H1,H2] = plotyy(date,MA(EmpGr(2:end),rolN,1)-mean(MA(EmpGr(2:end),rolN,1)),date,MktVol(1:end-1),'plot');
set(H1,'LineStyle','-','LineWidth',ms,'color','b')
set(H2,'LineStyle','-.','LineWidth',ms,'color','r')
set(AX(1),'XLim',[time0 time1])
set(AX(2),'XLim',[time0 time1])
xmin = -.02;
xmax = 0.02; %round(max(MA(IKq(2:end),rolN,1)),2);
ymin = -.5;
ymax = 0.5;% round(max(MktVol(1:end-1)),2);
set(AX(1),'YLim',[xmin xmax])
set(AX(2),'YLim',[ymin ymax])
set(AX(1),'YTick',[xmin 0.5*(xmin+xmax) xmax]); 
set(AX(2),'YTick',[ymin 0.5*(ymin+ymax) ymax]); 
set(AX(1),'ycolor','b','Fontsize',fs) 
set(AX(2),'ycolor','r','Fontsize',fs)
title('Employment growth','FontSize',fs) 
grid on

subplot(3,2,3)
[AX,H1,H2] = plotyy(date,MA(LiquidATGr(2:end),rolN,1)-mean(MA(LiquidATGr(2:end),rolN,1)), date,MktVol(1:end-1),'plot');
set(H1,'LineStyle','-','LineWidth',ms,'color','b')
set(H2,'LineStyle','-.','LineWidth',ms,'color','r')
set(AX(1),'XLim',[time0 time1])
set(AX(2),'XLim',[time0 time1])
xmin = -.25;
xmax = 0.25; %round(max(MA(IKq(2:end),rolN,1)),2);
ymin = -.5;
ymax = 0.5;% round(max(MktVol(1:end-1)),2);
set(AX(1),'YLim',[xmin xmax])
set(AX(2),'YLim',[ymin ymax])
set(AX(1),'YTick',[xmin 0.5*(xmin+xmax) xmax]); 
set(AX(2),'YTick',[ymin 0.5*(ymin+ymax) ymax]); 
set(AX(1),'ycolor','b','Fontsize',fs) 
set(AX(2),'ycolor','r','Fontsize',fs)
title('Cash growth','FontSize',fs) 
grid on

subplot(3,2,4)
[AX,H1,H2] = plotyy(date,MA(DivGr(2:end),rolN,1),date,MktVol(1:end-1) ,'plot');
set(H1,'LineStyle','-','LineWidth',ms,'color','b')
set(H2,'LineStyle','-.','LineWidth',ms,'color','r')
set(AX(1),'XLim',[time0 time1])
set(AX(2),'XLim',[time0 time1])
xmin = -.1;
xmax = 0.1; %round(max(MA(IKq(2:end),rolN,1)),2);
ymin = -.5;
ymax = 0.5;% round(max(MktVol(1:end-1)),2);
set(AX(1),'YLim',[xmin xmax])
set(AX(2),'YLim',[ymin ymax])
set(AX(1),'YTick',[xmin 0.5*(xmin+xmax) xmax]); 
set(AX(2),'YTick',[ymin 0.5*(ymin+ymax) ymax]); 
set(AX(1),'ycolor','b','Fontsize',fs) 
set(AX(2),'ycolor','r','Fontsize',fs)
title('Dividend growth','FontSize',fs) 
grid on


subplot(3,2,5)
[AX,H1,H2] = plotyy(date,MA(DebtGr(2:end),rolN,1)-mean(MA(DebtGr(2:end),rolN,1)),date,MktVol(1:end-1),'plot');
set(H1,'LineStyle','-','LineWidth',ms,'color','b')
set(H2,'LineStyle','-.','LineWidth',ms,'color','r')
set(AX(1),'XLim',[time0 time1])
set(AX(2),'XLim',[time0 time1])
xmin = -.05;
xmax = 0.05; %round(max(MA(IKq(2:end),rolN,1)),2);
ymin = -.5;
ymax = 0.5;% round(max(MktVol(1:end-1)),2);
set(AX(1),'YLim',[xmin xmax])
set(AX(2),'YLim',[ymin ymax])
set(AX(1),'YTick',[xmin 0.5*(xmin+xmax) xmax]); 
set(AX(2),'YTick',[ymin 0.5*(ymin+ymax) ymax]); 
set(AX(1),'ycolor','b','Fontsize',fs) 
set(AX(2),'ycolor','r','Fontsize',fs)
title('Debt growth','FontSize',fs) 
grid on

orient landscape

end
